
%Script to read crossover errors of individual trackline surveys, calculate
%error statistics
clear;
localdirectory = '/Users/manojnair/Downloads/';
%get data file names
%s = dir ('/Users/manojnair/projects/trackline_error/AllTotalField/');
%s = dir ('/Users/manojnair/projects/trackline_error/AllResidual/');
s = dir([localdirectory '*.cross']);

%# lon	lat	t_1	t_2	dist_1	dist_2	head_1	head_2	vel_1	vel_2	drt_1	drt_2	tz_1	tz_2	ptc_1	ptc_2	twt_1	twt_2	depth_1	depth_2	bcc_1	bcc_2	btc_1	btc_2	mtf1_1	mtf1_2	mtf2_1	mtf2_2	mag_1	mag_2	msens_1	msens_2	diur_1	diur_2	msd_1	msd_2	gobs_1	gobs_2	eot_1	eot_2	faa_1	faa_2	nqc_1	nqc_2

%loop through all surveys

%for i = 4:length(s),
for i = 1:length(s),
    
    
%A = textread (['/Users/manojnair/projects/trackline_error/AllTotalField/' s(i).name], '%f','headerlines', 4);
%A = textread (['/Users/manojnair/projects/trackline_error/AllResidual/' s(i).name], '%f','headerlines', 4);
[lon	lat	t_1	t_2,	dist_1	dist_2	head_1	head_2	vel_1	vel_2	drt_1...
    drt_2	tz_1	tz_2	ptc_1	ptc_2	twt_1	twt_2	depth_1	depth_2...
    bcc_1	bcc_2	btc_1	btc_2	mtf1_1	mtf1_2	mtf2_1	mtf2_2	mag_1...
    mag_2	msens_1	msens_2	diur_1	diur_2	msd_1	msd_2	gobs_1	gobs_2...
    eot_1	eot_2	faa_1	faa_2	nqc_1	nqc_2] =  textread([localdirectory s(i).name],'%f%f%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f','headerlines', 4);



B = reshape(A,[12,length(A)/12])';

s(i).rms = rms(B(:,11));
s(i).mean = mean(B(:,11));
s(i).median = median(B(:,11));
s(i).std = std(B(:,11));
s(i).meanlat = mean(B(:,2));
s(i).meanlong = mean(B(:,1));
s(i).medianlat = median(B(:,2));
s(i).medianlong = median(B(:,1));
s(i).ncorss = size(B,1);

% find robust statistics


[p,stats] = robustfit(1:length(B), B(:,11));
L = stats.w == 0;
s(i).robmean = mean(B(~L,11));
s(i).robrms = rms(B(~L,11));
s(i).robrejection = sum(L);
end;

%% write the data to a CSV file
fid = fopen('/Users/manojnair/projects/trackline_error/Total_Field_Error_2.csv','wt');
    fprintf(fid,'Survey_ID,Mean_Lat,Lean_Long,Median_Lat,Median_Long,N_Crosses,RMS,MEAN,STD,ROB_MEAN,ROB_MEDIAN,ROB_RMS,N_REJECTION,ABS_LAT,SQRT_ROB_RMS\n');

for i = 4:length(s),

    fprintf(fid,'%s,%f,%f,%f,%f,%d,%f,%f,%f,%f,%f,%f,%d,%f,%f\n', ...
    s(i).name, ...
    s(i).meanlat, ...
    s(i).meanlong, ...
    s(i).medianlat, ...
    s(i).medianlong, ...
    s(i).ncorss, ...
    s(i).rms, ...
    s(i).mean, ...
    s(i).median, ...
    s(i).std, ...
    s(i).robmean, ...
    s(i).robrms, ...
    s(i).robrejection, ...
    abs(s(i).medianlat),...
    sqrt(s(i).robrms));

end;
fclose(fid);


%% Read the crossover files with CM4 values.
%Script to read crossover errors of individual trackline surveys, calculate
%error statistics

clear;
close all;

% worldmap('world')
% load coast
% plotm(lat,long,'b','LineWidth',2);
% hold on;

% get the survey data from the xcel files
%[a,b] = xlsread('/Users/manojnair/projects/trackline_error/TotalFieldErrorGeomagLat_edited.xlsx');
[a,b] = xlsread('/Users/manojnair/projects/trackline_error/TotalFieldErrorGeomagLat.xlsx');
%localdirectory = '/Users/manojnair/projects/trackline_error/cm4_removed/crossoverErrors1/';
localdirectory = '/Users/manojnair/projects/trackline_error/cm4_removed/crossoverErrorsnew/';

% load Ap data for disturbance check
load /Users/manojnair/projects/trackline_error/Ap_1950_2012 Ap fday_Ap


%s = dir([localdirectory '*.cross']);
s = dir([localdirectory '*.mat']);




%loop through all surveys
for i = 1:length(s),
%for i = 29:29, % kr0112.mat is removed the data looks wiered
 
    
load([localdirectory s(i).name],'-mat');

% convert time to matlab
survey_fday_1 = datenum(t_1,'yyyy-mm-ddTHH:MM:SS');
survey_fday_2 = datenum(t_2,'yyyy-mm-ddTHH:MM:SS');


% get the mean Ap during the survey
s(i).mean_ap_1 = mean(interp1(fday_Ap, Ap, survey_fday_1));
s(i).mean_ap_2 = mean(interp1(fday_Ap, Ap, survey_fday_2));
s(i).mean_t1 = mean(survey_fday_1);
s(i).mean_t2 = mean(survey_fday_2);

s(i).rms = rms(crossErr);
s(i).mean = mean(crossErr);
s(i).median = median(crossErr);
s(i).std = std(crossErr);
[s(i).meanlat,s(i).meanlong] = meanm(lat,lon);
s(i).ncross = length(lon);
s(i).depth = median(depth_1);
s(i).diurnal = median(diur_1);


% arb_x_axis = [1:length(crossErr)];
% plot(arb_x_axis, crossErr, 'b.');
% hold on;

% find robust statistics

if length(crossErr) > 2,
% L = isnan(crossErr);
% lon(L) = [];
% crossErr(L) = [];

[p,stats] = robustfit(1:length(lon), crossErr);
L = stats.w == 0;
s(i).rms_without_cm4 = rms(mtf1_2(~L) - mtf1_1(~L));
s(i).robmean = mean(crossErr(~L));
s(i).robrms = rms(crossErr(~L));
s(i).robrejection = sum(L);
% plot(arb_x_axis(L), crossErr(L), 'r.');

else

s(i).robmean = NaN;
s(i).robrms = NaN;
s(i).robrejection = NaN;
end;





% get the magnetic latitude
path(path,'/Users/manojnair/projects/GIC/SHA_obs');
[theta_gm, phi_gm] = gg2gm_new(90 - s(i).meanlat,	s(i).meanlong);
s(i).maglat = - (theta_gm - 90);
% get the survey metadata
k = find(strcmp(b(2:end,1),[s(i).name(12:end-4) '_tCross'])==1);

if ~isempty(k)

s(i).survey_year = a(k,1);
s(i).gps = a(k,3);
s(i).residua= a(k,2);
s(i).oldrms = a(k,14);
else
s(i).maglat = NaN;
s(i).survey_year = NaN;
s(i).gps = NaN;
s(i).residua= NaN;
s(i).oldrms = NaN;
end;

%title([s(i).name(12:end-4) sprintf(', Mag Lat %4.2f, robrms = %4.2f, rms = %4.2f, rms_without_cm4= %4.2f',s(i).maglat,s(i).robrms,s(i).rms,s(i).rms_without_cm4)]);
%plotm(lat,lon,'r','LineWidth',2);
%plotm(s(i).meanlat, s(i).meanlong,'g*');
%textm(s(i).meanlat, s(i).meanlong,s(i).name(12:end-4),'FontSize',16);
% pause;
% cla;

end;



%%
data = load( '/Users/manojnair/projects/trackline_error/survey_apex_lat_long.txt');

plot(abs(data(:,4)),[s.rms_without_cm4],'r.','MarkerSize',30);
hold on;
plot(abs(data(:,4)),[s.robrms],'b.','MarkerSize',30);

set(gca,'FontSize',20)
xlabel('Magnetic latitude (ABS)');
ylabel('RMS (nT)');
legend('Original','CM4 Subtracted');
title('Crossover error as a function of magnetic latitude')

